El dengue es una enfermedad transmitida por la picadura de un mosquito que ocurre en lugares tropicales y subtropicales del mundo. Debido a que el vector de transmisión es un mosquito, la dinámica de la enfermedad se relaciona con variables climáticas como la lluvia y la temperatura. Aunque la relación con el clima es compleja, un gran número de científicos pronostica que el cambio climático puede aumentar los territorios en dónde es endémica la enfermedad lo cual tiene implicaciones importantes en materia de salud pública a nivel global.
De acuerdo a la Organización Mundial de la Salud (OMS):
“(…) Antes de 1970, solo nueve países habían sufrido epidemias de dengue grave. Ahora,la enfermedad es endémica en más de 100 países de las regiones de la OMS de África, las Américas, el Mediterráneo Oriental, Asia Sudoriental y el Pacífico Occidental. Las regiones más gravemente afectadas son las Américas, Asia Sudoriental y el Pacífico Occidental; de hecho en Asia se concentra aproximadamente el 70% de la carga mundial de la enfermedad.
Además de que el número de casos aumenta a medida que la enfermedad se propaga a nuevas zonas, incluida Europa, se producen brotes epidémicos de carácter explosivo. Europa ya se enfrenta a la posibilidad de brotes de dengue; la transmisión local se notificó por vez primera en Francia y Croacia en 2010 y se detectaron casos importados en otros tres países europeos. En 2012, un brote de dengue en el archipiélago de Madeira (Portugal) ocasionó más 2000 casos y se detectaron casos importados en otros 10 países europeos, además de Portugal continental. Actualmente se observan casos autóctonos casi anualmente en muchos países europeos. Entre los viajeros que regresan de países de ingresos bajos y medianos, el dengue constituye la segunda causa de fiebre más diagnosticada tras el paludismo.(…)
En 2019 se registró el mayor número casos de dengue jamás notificado en todo el mundo. Todas las regiones se vieron afectadas (…)” 1
Una mejor comprensión de la relación entre el clima y el dengue puede ayudar a las iniciativas de investigación en contra de la enfermedad y a una mejor distribución de los recursos económicos para la lucha contra de la propagación del vector.
La plataforma DrivenData 2 que tiene como objetivo incentivar la colaboración del público en proyectos de Data Science con impacto social se ha hecho eco de la iniciativa del gobierno norteamericano para luchar contra el dengue.
Los datos para esta competencia que patrocina la iniciativa “Predict the Next Pandemic Initiative”https://dengueforecasting.noaa.gov/ vienen de múltiples fuentes:
Los datos se anexarán a la presente investigación y se pueden conseguir en: https://www.drivendata.org/competitions/44/dengai-predicting-disease-spread/data/
Usando la data ambiental el objetivo del presente trabajo es predecir el número de casos totales de dengue reportados semanalmente en la ciudad de San Juan (Puerto Rico) y en la ciudad de Iquitos (Perú).
A continuación la lista de variables como las encontraremos en en el fichero:
City and date indicators
*city – City abbreviations: sj for San Juan and iq for Iquitos
*week_start_date – Date given in yyyy-mm-dd format
NOAA’s GHCN daily climate data weather station measurements * station_max_temp_c – Maximum temperature * station_min_temp_c – Minimum temperature * station_avg_temp_c – Average temperature * station_precip_mm – Total precipitation * station_diur_temp_rng_c – Diurnal temperature range
PERSIANN satellite precipitation measurements (0.25x0.25 degree scale) * precipitation_amt_mm – Total precipitation
NOAA’s NCEP Climate Forecast System Reanalysis measurements (0.5x0.5 degree scale)
Satellite vegetation - Normalized difference vegetation index (NDVI) - NOAA’s CDR Normalized Difference Vegetation Index (0.5x0.5 degree scale) measurements
Se cargan las librerías y se importa el conjunto de datos.
A continuación cargaremos los datos y los dividiremos de acuerdo a la ciudad, ya que suponemos que cada ciudad tiene un comportamiento distinto en cuanto a los casos de dengue y a sus condiciones climáticas. Es cierto que aunque ambas ciudades tienen un clima tropical, Iquitos es continental y San Juan se encuentra en una isla del Caribe.
##
## iq sj
## 35.71 64.29
##
## San Juan
## features: 24 entries: 936 labels: 936
##
## Iquitos
## features: 24 entries: 520 labels: 520
Un simple vistazo a la data nos revela que hay 24 columnas. De total del dataset un 36% corresponde a información sobre la ciudad de Iquitos y una mayoría de 64% a data sobre la ciudad de San Juan. Por lo pronto, como la variable objetivo es númerica, decimos que este será un modelo supervisado de regresión.
Una primera exploración de nuestros datos es necesaria para buscar patrones. Podemos observar el tipo de datos, el número de datos faltantes, las distribuciones de las variables, su correlación, y este análisis previo del dataset nos permitirá obtener una idea de las transformaciones que serán necesarias para su modelado, así como hacer hipótesis del tipo de modelo que podría ser útil. Haremos la exploración por separado para las dos ciudades apoyándonos en la visualización de los datos.
Empezamos con la ciudad de San Juan. A continuación algunos gráficos:
Conclusiones para ambas ciudades
De acuerdo a lo observado podemos extraer algunas conclusiones:
Tenemos 24 variables de las cuales 20 son númericas, 2 son enteros y 2 son de tipo carácter.
Las variables que se refieren a la vegetación tienen mayor proporción de missings. El caso más importante es la variable nvdi_ne con un 20% de datos faltantes para la ciudad de San Juan.
A simple vista, no se observa nada extraño en las distribuciones de los datos. No parecen haber outliers importantes.
Las variables ambientales están fuertemente correlacionadas entre ellas ya que se refieren en mucho casos al mismo fenómeno (lluvia, por ejemplo). Esto podría generar problemas pero de momento se trabajará con todas las variables.
El feature engineering es un proceso previo al modelado en donde generamos ajustes a las variables existentes o creamos variables nuevas a partir de las que nos fueron dadas con el objetivo de obtener mejores resultados en el modelo final de machine learning. A continuación la optimización del dataset.
En esta fase lo que vamos a hacer es crear algunas variables nuevas que ya vemos que pueden ser interesantes, por ejemplo el mes del año. Ya que es sabido que los fénomenos climáticos tienen cierto componente estacional.Así que puede haber una relación entre el mes, las condiciones atmosféricas y la proliferación de mosquitos.
testdengue<- as.data.table(test)
sj_train1 <- as.data.table(sj_train)
iq_train1 <- as.data.table(iq_train)
datfeat<- rbind(sj_train1, iq_train1)
datfeat<- rbind(datfeat, testdengue)Otra cosa que nos interesa es saber la distribución de las categorías para determinar cuáles son candidatas a one hot y cuáles a reemplazar por frecuencias.
unique(datfeat[ , .(.N), by = .(station_max_temp_c)][order(-N)])
unique(datfeat[ , .(.N), by = .(station_min_temp_c)][order(-N)])
unique(datfeat[ , .(.N), by = .(station_avg_temp_c)][order(-N)])
unique(datfeat[ , .(.N), by = .(station_precip_mm )][order(-N)])
unique(datfeat[ , .(.N), by = .(station_diur_temp_rng_c )][order(-N)])
unique(datfeat[ , .(.N), by = .(precipitation_amt_mm)][order(-N)])
unique(datfeat[ , .(.N), by = .( reanalysis_sat_precip_amt_mm )][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_dew_point_temp_k)][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_air_temp_k)][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_relative_humidity_percent )][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_specific_humidity_g_per_kg )][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_precip_amt_kg_per_m2)][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_max_air_temp_k )][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_min_air_temp_k)][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_avg_temp_k)][order(-N)])
unique(datfeat[ , .(.N), by = .(reanalysis_tdtr_k )][order(-N)])
unique(datfeat[ , .(.N), by = .(ndvi_se)][order(-N)])
unique(datfeat[ , .(.N), by = .(ndvi_sw )][order(-N)])
unique(datfeat[ , .(.N), by = .(ndvi_ne)][order(-N)])
unique(datfeat[ , .(.N), by = .(ndvi_nw )][order(-N)])Después de análisis:
1. Variables candidatas a transformación por frecuencias:
2. Variables que será necesario transformar por el último valor: station_avg_temp_c, station_diur_temp_rng_c , reanalysis_dew_point_temp_k, reanalysis_air_temp_k, reanalysis_relative_humidity_percent, reanalysis_specific_humidity_g_per_kg, reanalysis_precip_amt_kg_per_m2, reanalysis_avg_temp_k, ndvi_se, ndvi_sw, ndvi_ne, ndvi_nw
A continuación eliminamos variables que corresponden e imputamos los NAs por la moda en las variables. Es importante eliminar los NAS o imputarlos ya que los modelos de predicción generan error si hubiera datos faltantes.
datMod <- copy(datfeat)
#----Eliminamos variables que inspeccionamos y son redundantes.
datMod$week_start_date <- NULL
#-----------------------------------------------------
#names(datMod)
datMod_cl<- datMod
#Imputamos los NA que vimos en EDA. (por la moda)
getmode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
datMod_cl[is.na(datMod_cl$station_max_temp_c), "station_max_temp_c"] <- getmode(datMod_cl$station_max_temp_c)
datMod_cl[is.na(datMod_cl$station_min_temp_c), "station_min_temp_c"] <- getmode(datMod_cl$station_min_temp_c)
datMod_cl[is.na(datMod_cl$station_precip_mm ), "station_precip_mm"] <- getmode(datMod_cl$station_precip_mm )
datMod_cl[is.na(datMod_cl$precipitation_amt_mm ), "precipitation_amt_mm"] <- getmode(datMod_cl$precipitation_amt_mm )
datMod_cl[is.na(datMod_cl$reanalysis_sat_precip_amt_mm ), "reanalysis_sat_precip_amt_mm"] <- getmode(datMod_cl$reanalysis_sat_precip_amt_mm )
datMod_cl[is.na(datMod_cl$reanalysis_max_air_temp_k ), "reanalysis_max_air_temp_k"] <- getmode(datMod_cl$reanalysis_max_air_temp_k )
datMod_cl[is.na(datMod_cl$reanalysis_min_air_temp_k ), "reanalysis_min_air_temp_k"] <- getmode(datMod_cl$reanalysis_min_air_temp_k )
datMod_cl[is.na(datMod_cl$reanalysis_tdtr_k ), "reanalysis_tdtr_k"] <- getmode(datMod_cl$reanalysis_tdtr_k )#Variables que será necesario transformar por el último valor:** station_avg_temp_c, station_diur_temp_rng_c , reanalysis_dew_point_temp_k, reanalysis_air_temp_k, reanalysis_relative_humidity_percent, reanalysis_specific_humidity_g_per_kg, reanalysis_precip_amt_kg_per_m2, reanalysis_avg_temp_k, ndvi_se, ndvi_sw, ndvi_ne, ndvi_nw
datMod_cl$ndvi_ne %<>% na.locf(fromLast = TRUE)
datMod_cl$ndvi_sw %<>% na.locf(fromLast = TRUE)
datMod_cl$ndvi_se %<>% na.locf(fromLast = TRUE)
datMod_cl$ndvi_nw %<>% na.locf(fromLast = TRUE)
datMod_cl$station_avg_temp_c %<>% na.locf(fromLast = TRUE)
datMod_cl$station_diur_temp_rng_c %<>% na.locf(fromLast = TRUE)
datMod_cl$reanalysis_dew_point_temp_k %<>% na.locf(fromLast = TRUE)
datMod_cl$reanalysis_air_temp_k %<>% na.locf(fromLast = TRUE)
datMod_cl$reanalysis_relative_humidity_percent %<>% na.locf(fromLast = TRUE)
datMod_cl$reanalysis_specific_humidity_g_per_kg%<>% na.locf(fromLast = TRUE)
datMod_cl$reanalysis_precip_amt_kg_per_m2 %<>% na.locf(fromLast = TRUE)
datMod_cl$reanalysis_avg_temp_k%<>% na.locf(fromLast = TRUE)Las variables que anteriormente se identificaron como candidatas a transformar por frecuencias se modifican de catégoricas a numéricas de esta manera y se eliminan las originales.
# Voy a convertir a Frequencia las variables
datMod_cl[ , fe_reanalysis_sat_precip_amt_mm := .N , by = .(reanalysis_sat_precip_amt_mm )]
datMod_cl[ , fe_reanalysis_dew_point_temp_k := .N , by = .(reanalysis_dew_point_temp_k )]
datMod_cl[ , fe_reanalysis_air_temp_k := .N , by = .(reanalysis_air_temp_k )]
datMod_cl[ , fe_reanalysis_relative_humidity_percent := .N , by = .(reanalysis_relative_humidity_percent )]
datMod_cl[ , fe_reanalysis_specific_humidity_g_per_kg := .N , by = .(reanalysis_specific_humidity_g_per_kg)]
datMod_cl[ , fe_reanalysis_precip_amt_kg_per_m2 := .N , by = .(reanalysis_precip_amt_kg_per_m2 )]
datMod_cl[ , fe_reanalysis_max_air_temp_k := .N , by = .(reanalysis_max_air_temp_k )]
datMod_cl[ , fe_reanalysis_min_air_temp_k := .N , by = .(reanalysis_min_air_temp_k)]
datMod_cl[ , fe_reanalysis_avg_temp_k := .N , by = .(reanalysis_avg_temp_k )]
datMod_cl[ , fe_reanalysis_tdtr_k := .N , by = .(reanalysis_tdtr_k )]
datMod_cl[ , fe_ndvi_se := .N , by = .(ndvi_se)]
datMod_cl[ , fe_ndvi_sw := .N , by = .(ndvi_sw )]
datMod_cl[ , fe_ndvi_ne:= .N , by = .(ndvi_ne)]
datMod_cl[ , fe_ndvi_nw := .N , by = .(ndvi_nw )]
datMod_cl[ , fe_precipitation_amt_mm := .N , by = .(precipitation_amt_mm )]
datMod_cl[ , fe_station_diur_temp_rng_c := .N , by = .( station_diur_temp_rng_c )]
datMod_cl[ , fe_station_avg_temp_c := .N , by = .(station_avg_temp_c )]
datMod_cl[ , fe_station_precip_mm := .N , by = .(station_precip_mm)]
datMod_cl[ , fe_station_max_temp_c := .N , by = .(station_max_temp_c)]
datMod_cl[ , fe_station_min_temp_c:= .N , by = .(station_min_temp_c)]
to_rem <- c('station_avg_temp_c', 'station_precip_mm', 'station_diur_temp_rng_c','precipitation_amt_mm', 'reanalysis_sat_precip_amt_mm', 'reanalysis_dew_point_temp_k', 'reanalysis_air_temp_k','reanalysis_relative_humidity_percent', 'reanalysis_specific_humidity_g_per_kg', 'reanalysis_precip_amt_kg_per_m2','reanalysis_max_air_temp_k', 'reanalysis_min_air_temp_k' , 'reanalysis_avg_temp_k', 'reanalysis_tdtr_k', 'ndvi_se','ndvi_sw', 'ndvi_ne', 'ndvi_nw','station_max_temp_c', 'station_min_temp_c' )
datMod_cl[ , (to_rem) := NULL]
#datMod_cl %>%
# select(-status_group)%>%
#glimpse()Ahora, vuelvo a separar el data set a su forma original pero con las transformaciones hechas en ambos. Además cambio el nombre de los niveles de la variable objetivo para que sean nombres “aceptables” por R y no den error en el fit.
Por último, compruebo que se mantienen aproximadamente los valores en porcentaje de cada categoría de la variable target.
#Divido de nuevo los dataset en train y test
dat_train<- split(datMod_cl1,findInterval(1:nrow(datMod_cl1), 1457))
df_train<- dat_train[[1]]
df_test<- dat_train[[2]]
df_train<- merge(df_train,train_labels, by=c("city", "year", "weekofyear"))
df_train$city <- as.factor(df_train$city)
levels(df_train$city)
datMod_cl1<-df_train
# Compruebo de nuevo proporciones del target
round(prop.table(table(datMod_cl1$city))*100, 2)Una vez procesados los datos, divido mi conjunto de entrenamiento en test y train. En este punto fijo una semilla para que sean comparables los resultados de los fit.
seed<-7
set.seed(seed)
validationIndex <- createDataPartition(datMod_cl1$total_cases, p = 0.80, list = FALSE)
# select 20% of the data for validation
my_test <- datMod_cl1[-validationIndex,]
# use the remaining 80% of data to training
my_train <- datMod_cl1[validationIndex,]
iq_train2 = my_train %>% filter(city == 'iq')
sj_train2 = my_train %>% filter(city == 'sj')
iq_test2 = my_test %>% filter(city == 'iq')
sj_test2 = my_test %>% filter(city == 'sj')
library(janitor)
my_train<- clean_names(my_train)
iq_train2 <- clean_names(iq_train2)
sj_train2 <- clean_names(sj_train2)
iq_test2 <- clean_names(iq_test2)
sj_test2 <- clean_names(sj_test2)
my_test<- clean_names(my_test)A continuación utilizo validación cruzada.
Una vez tenemos los conjuntos de entrenamiento y de test (creados por nosotros) vamos al modelo. Haremos un modelo de Random Forest. Un modelo Random Forest está formado por un conjunto de árboles de decisión individuales, cada uno entrenado con una muestra ligeramente distinta de los datos de entrenamiento generada mediante bootstrapping). La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.
library(MLmetrics)
library(ranger)
set.seed(seed)
fit <- ranger(
total_cases ~. ,
data = sj_train2,
num.trees = 100,
importance = 'impurity',
write.forest = TRUE,
min.node.size = 1,
verbose = TRUE,
classification = FALSE,
)Veo resultados del modelo
## Ranger result
##
## Call:
## ranger(total_cases ~ ., data = sj_train2, num.trees = 100, importance = "impurity", write.forest = TRUE, min.node.size = 1, verbose = TRUE, classification = FALSE, )
##
## Type: Regression
## Number of trees: 100
## Sample size: 746
## Number of independent variables: 24
## Mtry: 4
## Target node size: 1
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 922.5095
## R squared (OOB): 0.6984497
## Length Class Mode
## predictions 746 -none- numeric
## num.trees 1 -none- numeric
## num.independent.variables 1 -none- numeric
## mtry 1 -none- numeric
## min.node.size 1 -none- numeric
## variable.importance 24 -none- numeric
## prediction.error 1 -none- numeric
## forest 8 ranger.forest list
## splitrule 1 -none- character
## treetype 1 -none- character
## r.squared 1 -none- numeric
## call 10 -none- call
## importance.mode 1 -none- character
## num.samples 1 -none- numeric
## replace 1 -none- logical
## dependent.variable.name 1 -none- character
Puedo incluso ver las variables importantes
Creo el gráfico de importancia de las variables.
Intento mejorar el resultado de mi modelo ranger cambiando los hiperparámetros.
library(MLmetrics)
library(ranger)
set.seed(seed)
# Barrido con diferente número de arboles.
val_trees <- c(100, 150, 200, 250)
for (i in val_trees) {
print(i)
fit_rg <- ranger(
total_cases~. ,
data = sj_train2,
num.trees = i,
importance = 'impurity',
write.forest = TRUE,
min.node.size = 1,
verbose = TRUE,
classification = FALSE
)
valor_pred <- predict(fit_rg, data = sj_test2)
fit_mae <- MAE(y_pred = valor_pred$predictions, y_true = sj_test2$total_cases)
print(fit_mae)
}## [1] 100
## [1] 15.67399
## [1] 150
## [1] 15.71255
## [1] 200
## [1] 15.22339
## [1] 250
## [1] 15.92086
El mejor valor de num.trees es de 100 con el que se consigue un MAE de 15.64.
library(MLmetrics)
library(ranger)
set.seed(seed)
fit_iq <- ranger(
total_cases ~. ,
data = iq_train2,
num.trees = 100,
importance = 'impurity',
write.forest = TRUE,
min.node.size = 1,
verbose = TRUE,
classification = FALSE,
)Veo resultados del modelo
## Ranger result
##
## Call:
## ranger(total_cases ~ ., data = iq_train2, num.trees = 100, importance = "impurity", write.forest = TRUE, min.node.size = 1, verbose = TRUE, classification = FALSE, )
##
## Type: Regression
## Number of trees: 100
## Sample size: 422
## Number of independent variables: 24
## Mtry: 4
## Target node size: 1
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 102.658
## R squared (OOB): 0.1869907
## Length Class Mode
## predictions 422 -none- numeric
## num.trees 1 -none- numeric
## num.independent.variables 1 -none- numeric
## mtry 1 -none- numeric
## min.node.size 1 -none- numeric
## variable.importance 24 -none- numeric
## prediction.error 1 -none- numeric
## forest 8 ranger.forest list
## splitrule 1 -none- character
## treetype 1 -none- character
## r.squared 1 -none- numeric
## call 10 -none- call
## importance.mode 1 -none- character
## num.samples 1 -none- numeric
## replace 1 -none- logical
## dependent.variable.name 1 -none- character
Puedo incluso ver las variables importantes
Creo el gráfico de importancia de las variables.
#### - Evalúo modelo ranger
valor_pred2 <- predict(fit_iq, data = iq_test2)
#table(my_test$total_cases, valor_pred$predictions)
# Usando el paquete MLmetrics puedo calcular diferentes métricas
library(MLmetrics)
MAE(y_pred = valor_pred2$predictions, y_true = iq_test2$total_cases)## [1] 4.825011
Intento mejorar el resultado de mi modelo ranger cambiando los hiperparámetros.
library(MLmetrics)
library(ranger)
set.seed(seed)
# Barrido con diferente número de arboles.
val_trees <- c(100, 150, 200, 250)
for (i in val_trees) {
print(i)
fit_rg_iq <- ranger(
total_cases~. ,
data = iq_train2,
num.trees = i,
importance = 'impurity',
write.forest = TRUE,
min.node.size = 1,
verbose = TRUE,
classification = FALSE
)
valor_pred3 <- predict(fit_rg_iq, data = iq_test2)
fit_mae <- MAE(y_pred = valor_pred3$predictions, y_true = iq_test2$total_cases)
print(fit_mae)
}## [1] 100
## [1] 4.825011
## [1] 150
## [1] 4.583321
## [1] 200
## [1] 4.773587
## [1] 250
## [1] 4.714154
El mejor valor de num.trees es de 250 con el que se consigue un MAE de 4.64.
Utilizo la data de test de la plataforma para obtener predicciones que subir al concurso. Para ello empleo el modelo ranger con grid search.
Vemos los resultados de la plataforma:
Los resultados no son tan buenos como esperabamos de acuerdo al entrenamiento. Sugerimos a futuras investigaciones sobre este problema tratar con distintos modelos o con una selección de variables distinta.